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1. Introduction 

Networks of coupled oscillators are used to describe a variety of systems in science 
and engineering, such as Josephson junction arrays, generators in power plants, firefly 
populations, and heart pacemaker cells. Of particular interest are solutions in which the 
network, or subpopulations within the network, oscillate synchronously. The analysis 
of the stability of and transition to a synchronous state can be very complex and has 
received much attention [1, 2, 3]. 

Recent applications to nanoelectromechanical systems (NEMS) [4], and beam 
steering devices in telecommunications [5], showed that important advances can be made 
by studying these problems perturbatively. It is therefore essential to have appropriate 
mathematical tools for such an analysis. We propose a perturbative method, based on 
normal form techniques [6, 7, 8], which is in many respects superior to those commonly 
used to study synchrony in oscillator networks. 

The method is intuitive and helps us distinguish between contributions to the 
dynamics arising from the network configuration and the internal structure of individual 
oscillators. Because of this it is possible to carry out calculations without having 
to specify the nonlinearities explicitly. Also, the approach based on normal forms is 
rigorous, and the validity of approximations is known a priori. On the other hand, for 
methods commonly used in the physics literature, mathematical justification is often 
non-trivial, and must be performed a posteriori [8, 9]. 

There are a number of other advantages that the normal form approach brings to the 
table. The approximating equations to the original system are obtained by examining a 
collection of algebraic conditions. This procedure can be formulated in an algorithmic 
form and automated, which is of particular importance when approximations of higher 
order in the small parameter are needed. Computer codes for determining normal forms 
to any order already exist for some problems in celestial mechanics. 

Standard methods, such as averaging, usually require center manifold reduction to 
be performed first [10]. We will show that the center manifold reduction is obtained 
naturally in the normal form of the equations of motion. Moreover, the change of 
coordinates leading to the normal form can be used to approximate the center manifold, 
the invariant fibration over the center manifolds, and a number of nearly conserved 
quantities of the equations to any order. As a result, the normal form method offers a 
deeper insight into the geometric structure of the approximating equations. 

Related approaches can be found in the literature [11,12]. In particular, the method 
of normal forms has been used in [13] to obtain reduced equations for oscillators close 
to a bifurcation. Our approach differs in that we do not consider only small amplitude 
oscillations, but general weakly nonlinear oscillators. Moreover, the coupling in the 
present case results in negative eigenvalues in the linear part of the vector field. 

In this paper we look at systems of globally coupled identical oscillators illustrated 
in Figure 1. This configuration is commonly used for wave generators in order to increase 
the output power (See e.g. [] 1]). When oscillators are synchronized the power of emitted 
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Figure 1. Globally coupled oscillators. Load parameters are rescaled with respect to 
the number of oscillators, so the values of load resistance, capacitance and inductance 
are NR, C/N and NL, respectively. 

waves scales as a square of the number of coupled units. Therefore, it is important to 
determine couplings that lead to synchronous behavior. Our analysis results in a general 
expression for the onset of synchronization in the network, and we recover recent results 
for an array of van der Pol oscillators [15] as a special case. Furthermore, we find, 
somewhat surprisingly, that the coupling can induce synchronous oscillations even in 
a network of weakly nonlinear systems which are unstable, and do not oscillate when 
uncoupled. The method itself can be easily extended to treat more complex networks 
and other types of coherent solutions. In order to keep our presentation streamlined we 
do not discuss these problems here. 

This work is motivated by previous studies of synchrony in Josephson junction 
arrays [16, 17, 18, 19], where a series of junctions was shunted with an RLC load 
(Fig. 1). Dhamala and Wiesenfeld introduced a heuristic, perturbative method in 
which an approximate stroboscopic map was constructed [19]. For an appropriately 
chosen strobing period T the synchronous solution corresponds to a fixed point of this 
stroboscopic map, and the stability of the synchronous state can be determined from the 
eigenvalues of its linearization. A remarkable consequence of this approach is a unified 
synchronization law for capacitive and noncapacitive junctions, two cases which were 
believed to have different dynamics. An extension of the method, and an application to 
the study of synchrony in an array of van der Pol oscillators was given in [15]. 

The stroboscopic map approach, like most classical singular perturbation methods, 
consists of identifying and taming secular terms in the naive approximating solution of 
the weakly nonlinear system. The general structure of the reduced equation obtained 
using this approach is difficult to know before the calculations are carried out. On the 
other hand, the normal form method enables us to carry out calculations without having 
to specify the nonlinearity explicitly, or calculate the approximate strobing time T, and 
hence study a much broader class of problems. 

The paper is organized as follows: In Section 2 we briefly review the theory of 
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normal forms and illustrate it in the case of a van der Pol oscillator. These ideas are 
applied in Section 3 to compute the normal form of the equation describing a network of 
weakly nonlinear oscillators. The stability of the synchronous solution in this network 
is analyzed in Section 4. Several examples illustrating these ideas are given in Section 
5. Finally, in Section 6 we discuss possible extensions of our work to different networks. 



2. Normal Form Analysis 

In this section we give outline the application of the method of normal forms to the 
analysis weakly nonlinear systems in a general setting. The method has been discussed 
in [7, 6], and a mathematical analysis is given in [20, 8]. Consider a weakly nonlinear 
ODE of the form 

x' = Ax + e ^ c„,ix"ei = Ax + ef (x), (1) 

ct.i 

where x G M", A is an n x n constant matrix, and is the i-th unit vector. We use 
standard multi-index notation, so that a G and x" = We emphasize 

that, in contrast with local normal form theory, f(x) may contain any monomial in x, 
including linear terms. 

A goal of normal form theory is to remove terms of 0{e) in equation (1) by a 
near-identity change of variables 

x = y + eg(y), g:M"^M", (2) 
where g(y) is a polynomial. In terms of the new variables (1) becomes 



r' = Ay + e[Y, c«,.y"e. - [A, g](y) j + 0{e' 



(3) 



where the Lie bracket [A, g](y) equals Dgiy^Ay — y4g(y). To remove the nonlinear 
terms at Oie) in (3), we need to solve the equation 

[Ag](y) = 5^c«,,,y"e,. (4) 

Since this equation is linear in g, it is equivalent to the finite family of equations 

[^,fi'a,i](y) = c«,jy"ej. (5) 

If A is diagonal, the eigenvectors of [A, ■] are the homogeneous polynomials, since 
[A, y"e,] = A„,»y"ei, 

where 

Aa,j = ^ ttfcAfc - \i = {a, X) - \i, (6) 

k 

and Aj are eigenvalues of matrix A. It follows that, if A^^i ^ 0, equation (5) has the 
solution 

gocAY) = T — y e*. 
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On the other hand, if Aq, j = 0, equation (5) does not have a solution, and we say that 
the monomial y"ej is resonant. Therefore, only nonlinear monomials y"ej such that 
Aa,i 7^ can be removed at first order in e by a near identity coordinate change of the 
form (2). 

In particular, we can split the nonlinearity f(y) in (1) into a resonant part 
f^(y) = ZlA„,,=o^«:iy"^i ^ nonresonant part f^'^(y) = Z^A^^i^o ^a.^y^ej, so that 
f(y) = f'^(y) + f^'^(y). By setting g(y) = Ea„,,^o fi'"-i(y)' change of coordinates 
(2) leads to the equation 

y' = Ay + 6f'^(y) + 0(e2). (7) 

We emphasize that to obtain the normal form of equation (1) to 0{e), we simply 
identify and remove all resonant terms, that is all monomials comprising f'^'^(x). The 
preceding argument shows that this can be done at the expense of introducing terms 
of 0{e^) into the equation. If the O(e^) terms are neglected in (7), a truncated normal 
form is obtained. To continue this process and obtain normal forms to higher order in 
e, it is necessary to compute the 0{e^) terms that are introduced at this step explicitly. 
This is equivalent to the observation that the computation of a local normal form near 
a singular point to second order may affect the cubic terms, and the computation needs 
to be carried out order by order (see [9], for instance). 

The following theorem shows that the truncated normal form provides a good 
approximation to the original equations 

Theorem 1 /s'/ Consider the ordinary differential equation 
x' = Ax + e^/„,,(t)x"ei, x(0)=xo 

where A is a diagonal, and has eigenvalues with non-positive real part. Construct the 
truncated normal form 

y' = Ay + ef'^(y), 

and let x(0) = y(0). Then there is aT = T(x(0)) > such that the solutions of the two 
equations satisfy |x(t) — y(t)| = 0{e) for all t G [0,T/e] and e sufficiently small. 

Remark 1 When the truncated equations contain hyperbolic invariant structures, 
the conclusions of Theorem 1 often hold for all time along their stable directions 
[21]. Normal form theory can be extended to study nondiagonalizable matrices [9], 
nonhomogeneous equations with nonlinearities that are not finite sums of monomials, 
and higher order approximations in e [8, 22]. The main ideas presented here are similar 
in these cases. We present the simplest case in order to keep technicalities at a minimum. 

Note that, by construction [A, i^] = 0. It follows that the truncated normal 
form is equivariant under the flow of the unperturbed equation x' = Ax. Moreover, 
the resonant monomials and hence the structure of the truncated normal form, are 
completely determined by the eigenvalues of A. This allows us to prove the following 
Proposition which will be useful in the following sections. 
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Proposition 1 Suppose that matrix A in (1) is diagonalizahle, has m purely imaginary 
eigenvalues Aj, and that other eigenvalues Vi have negative real part. The system (1) 
then can be written as 

x'l = AiXi + ehi(xi,X2) 
X2 = A2X2 + eh2(xi, X2) 

where Ai = diag{Xi, . . . , Xm), A2 = (z^m+i, • • • ^n.)- If the nonlinear terms hi and h2 are 
polynomials, the truncated normal form of this system has the general form 

yi = Aiyi + eh^(yi) (8) 

y'2 = My 2 + (yi , y2) (9) 

where (yi, 0) = 0. 

Proof: A term xf x^e^ with i < m is resonant if Ac,,i = ^Jli Q^i-^j + Sj=m+i ^i^i ~ 
Xi = 0. Since the eigenvalues z/j have negative real part and enter the sum with the 
same sign, this condition can be satisfied only if all (3j = 0. 

Similarly, a term x"ej for 2 > m is resonant only if Aq, j = Yl^Li '^j^j ~ = 0. This 
equation cannot hold. Hence all resonant monomials in (9) contain a nonzero power of 
y2j for some j, and evaluate to when y2 = 0. □ 

Although Proposition 1 is simple to prove, it says much about the structure of the 
truncated normal form. The fact that h2 (yi, 0) = means that the hyperplane y2 = 
is invariant under the fiow of (8-9). In fact, the hyperplane y2 = is the center manifold 
of this system, and hence (8) trivially provides the reduction of the truncated normal 
form equation to the center manifold. Therefore, it is unnecessary to compute the center 
manifold explicitly to obtain the reduced equations. 

Furthermore, since y2 does not occur on the right hand side of (8), the fibration 
given by yi = const, is also invariant under the fiow. To obtain an 0{e) approximation 
of the center manifold, and the invariant fibration over the center manifold in the 
original coordinates, it is sufiicient to invert the near identity transformation (2) used in 
obtaining the normal form equation. In fact, there are typically other easily identifiable 
quantities that are conserved by the fiow of the truncated normal form, and provide 
adiabatic invariants for the original equations [9, Chapter 5]. 

2.1. Example: Van der Pol Oscillator 

As a simple, illustrative example, and to introduce results that will be used in section 
5, we consider the van der Pol equation 

x" - ex' {I - x^) + X = 0. (10) 

We can rewrite system (10) in the variables 

z = X + ix\ z = X — ix (11) 
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to obtain: 

z' = - iz + -{Az - Az - - z^z + zz^ + z^), (12) 
8 

z' = iz- ^{Az -Az-z^- z^z + zt + f"). (13) 
8 

This system is in the form that can be analyzed using the ideas discussed in Section 2. 
The eigenvalues of the linear part are Ai = —i and A2 = i- The resonant terms can be 
computed using the condition Aq, j = where Aq, j is defined in equation (6). 



Table 1. Resonance condition. 
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As noted in the discussion following the derivation of equation (7), the truncated 
normal can be obtained simply by removing the resonant monomials in (12-13). From 
Table 1 we find that the resonant terms in (12) are z and z^z, and the resonant terms 
in (13) are their complex conjugates z and expected. Therefore, there exists a 

near identity change of coordinate in which (12-13) have the form 

z' = -tz + -z- -z^z + C(e2), (14) 
2 8 

and it's complex conjugate. 

As noted in Section 2, equation (14) is equivariant under the flow of the unperturbed 
equation, which is a pure rotation of the real plane. It follows that the right hand side of 
any weak perturbation of the equation z' = —iz, z! = iz cannot depend on the angular 
variable when expressed in polar coordinates. Indeed, (14) takes the form 

R' = (^1 - \R^y e' = i. 

in polar coordinates. The same procedure can be used to obtain higher order normal 
forms, see [8]. 

Remark 2 Strictly speaking, new coordinates are introduced in obtaining equation (14)- 
To keep notation at a minimum, we name these new variables z and z as well. A similar 
convention is used in the rest of the paper. 

3. Globally coupled networks 

In this section we use the normal form method to study a network of identical, weakly 
nonlinear oscillators, described by the equation x'- +Xi + eh{xi, x'^) = when uncoupled. 
Here e is a small parameter, and the nonlinear term h{xi,x'j) is assumed to be a 
polynomial. The elements in the network are globally coupled by a linear load (Fig. 
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1). The coupling is weak and of the same order as the nonhnearity. While this is not 
the most general example of a globally coupled network, these assumptions have been 
chosen to simplify the presentation, and can be relaxed. Equations of motion for this 
system can be written as 

x'l + eh{xk,x'^,) + Xk = q' (15) 



N 



/iig + Ai2g + q = e 



(16) 



Here we follow the notation introduced in [18, 16, 19, 15], where fii and /i2 are control 
parameters, and can be understood as the inductance and the resistance of the coupling 
load, respectively. The goal of our calculation is to find load parameters that will ensure 
synchrony in the network. 

In order to bring the equations of motion to normal form, we first have to diagonalize 
their linear parts. To make the procedure more intuitive we introduce complex variables 
Zk = + ix'^ and Zk = Xk — ix'f^, and denote q' = p. The system (15-16) then becomes 



where f{z, z) 



z'k = -izk + ip- ief{zk, Zk) 
4= iZk - ip + ief{zk,Zk) 
= h[{z + z)/2,{z - z)/2], and 
p 



^2 1 , 

— p q + 



N 



This system can be written in matrix form as 

z' = Az + ief (z) 

where z = {zi, ^i, 2:2 ... , Zn, q,p), 
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(22) 



.£^2. / 
Ml ' 



and 
f(z) 



T 



-/(^l, Zx), f{Zi, Zi), -f{z2, Z2),..., f{zN, Zn), 0, 



1K> ^ — ^ , _ 



(23) 



Here denotes the transpose, and f(z) is a column vector. 
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Finally, we diagonalize the matrix A, by introducing another coordinate change 
w = B~^z, where B~^AB = Da is diagonal. In the new coordinates the equations of 
motion have the form 

w' = B-^Az + ieB-^iz) = Daw + ieB-^i{Bw), (24) 

The matrix B can be computed using elementary linear algebra. Its actual form is not 
of interest, and we therefore suppress it. In component form (24) becomes 

= - iwk - iefiwk + biu + b2V, Wk + hu + b2v) (25) 

+ e^e^' 5^K- + + ^(^2 + h) + uibi + h)] 
j 

w'^ = iWk + ief{wk + biu + b2V + Wk + biu + b2v) (26) 

+ ^^e"'' ZlK- + + v{b2 + h) + u{bi + 6i)] 
i 

u' = l'2N+l u (27) 

- i{l - J2^Wj + Wj + v{b2 + 62) + u{bi + 61)] 

j 

V' = V2N+2 V (28) 



- ^(1 + C)^^ 5^K- + ^3 + <b2 + 62) + U{bi + 61)] 



where 61,2 = 2/ii/ (2/ii + z/i2 =F ZA/7if^^4/i^), C = /i2/ ^/J^^^Mm., Z = a/(1 - /ii)^ + /ij 
is the impedance, and 6 = arcsin(/i2/2') is the phase shift on the load. Variables u and 
V are obtained from q and p. Note that matrix A has 2N purely imaginary eigenvalues 
Afc = —h ^k+i = iyk = 1, 3, ... , 2N — 1, corresponding to the first 2N entries on the 
diagonal of B^^AB. The last two eigenvalues 



-fl2 ± V/ij - 4/ii 

i^2N+i,2 = i-^yj 

2/Ui 

have negative real part. 

Since the linear part of (25-28) is diagonal, it is now straightforward to identify 
nonresonant terms, as discussed in Section 2. First, from Proposition 1, we find that all 
terms containing powers of m or f can be removed at 0{e) in the equations for Wk and 
Wk using a near identity change of coordinates, so (25-26) reduce to the equations on 
the center manifold 

w',, = -iWk - iefiwk, Wk) + ^^^e'^ J^^^, + wj) + 0{e^), (30) 

j 

K 



w'k= iWk + tef{wk:Wk) + e^^e '^Y.^w, + w,) + 0{^). (31) 

j 



While the load variables do not enter the final approximating equations, the form of 
the load equation was important in obtaining the diagonalization, and will therefore be 
reflected in the final approximation. 

We can use the normal form to compute the transients, as well as the asymptotic 
state of a solution. However, since we are interested in the stability of the synchronous 
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state, we only consider the reduced equations on the center manifold u = 0, v = 0, and 
therefore drop equations (27-28) from further consideration. 

The normal form of (25-26) is obtained by computing the remaining resonant 
terms. The ones due to the nonlinearity /, are determined as follows: Let be the 
resonant part of / in the equation for the individual oscillators, w' = —iw — ief{w,w), 
and /? the resonant part in w' = iw + ief{w,w). A simple computation shows that 
fw — '^k'P^iwkWk) and = Wk4>^{wkWk), where (p^ is a polynomial in ww. Generally, 
the coefficients of (p^ are complex. 

It remains to determine the resonant terms that are due to the coupling term 
Kc^^ / {2ZN)^-{wj + Wj). Obviously, monomials Wj are resonant in (30), while 
monomials Wj are resonant in equations (31). Keeping only the resonant terms in the 
(30-31), we obtain the normal form to 0{e): 



wi = -iwk + e ( Wk(t>''{wkWk) + ^e*' | + 0{e') (32) 



w'l. 



iwk + e[ Wk^''{wkWk) + i^e''' 5^ ) + O(e'). (33) 



2ZN 

j 



Since the normal form equations are equivariant under the flow of the unperturbed 
system, it is again natural to rewrite them in polar coordinates, Wk = rkCT^^^- We obtain 



evkRivk) + 6^ - + + ^(^') (34) 

i 

Q'^ = \- eQivk) - ^ ^ sin(^, - 6, +6) + 0{e'), (35) 

j 

where R = {(j)^ + cf)^) / 2 is the real, and = [cj)^ — cf)^) / 2i the imaginary part of cf)^. Note 
that the fact that the right hand side of (34) depends only on the phase differences is a 
consequence of the equivariance of the truncated normal form under the transformation 
e -^e + C, where 6 = {Oi, 6*2, ... , 9n) and C = (c, c, . . . , c). 

Remark 3 We could also use a center manifold reduction to remove the variables u and 
V from equations (25-26) [10]. This would add another step to the calculation. If the 
transient dynamics of initial states off the center manifold is of interest, in addition it 
is necessary to compute the stable fibration over the center manifold. The normal form 
method considerably simplifies these computations. As noted in Section 2, the truncated 
normal form is a skew product, and provides both the reduction of the equations to the 
center manifold, and an approximation for the flow in the transversal direction [23]. 



4. Existence and stability of the synchronous state 

The truncated normal form obtained from (34-35) by neglecting 0(e^) terms is much 
easier to analyze than the original equations (15-16). It is straightforward to flnd the 
inphase solution for the truncated system and its stability analytically. Furthermore, if 
the truncated system has a stable synchronous solution, so does the full system (15-16). 
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Due to all-to-all coupling, the truncated normal form given by (34-35) is equivariant 
under all permutations of the oscillators, as well as the the action of the group given 
hj 9i — s> 9i + C for all i. As a consequence, a number of phase locked states are 
forced to exist [12]. In this section we consider the stability of the synchronous state 
Ti = r2 = . . . risf = r, and 6i = 62 = . . . = 6^ ■ The stability of other phase locked states 
can be analyzed similarly. 

Assuming that = r, and 6k — 9j = for all k and j, then equation (34) implies 
that r must satisfy 

rR{r) + r—— cos 6 = 0. (36) 

Remark 4 The amplitudes oj the uncoupled oscillators are given by solutions of R{r) = 
0, while in a synchronously oscillating network the amplitudes of the oscillators are given 
as solutions of (36) in terms of the coupling strength k, and properties of the coupling 
load. It is possible that R{r) = has only r = as a solution, while (36) has nonzero 
solutions for certain values of k, Z, 6. In such examples, the uncoupled systems do not 
oscillate, while the network can exhibit synchronous oscillations (see Section 5.3). 

Due to the equivariance of the system, the Jacobian is constant along the 
synchronous solution. Therefore the Floquet exponents equal the eigenvalues of the 
Jacobian. At the synchronous solution determined by 9k — 6j = and (36) they equal 

Ai =0 (37) 

A2 = erR'{r) (38) 

A„,n+i = erR'{r) - e— cos b (39) 

Zj 

2 



±eJ{rR'{r)y ~ [^sind^ +2r&{r)^sm6 



where n = 3, 5, 7, . . . , 2N — 1. If r satisfies (36) and the Aj>i have non-zero real part 
then the system obtained by truncating (34-35) has a weakly hyperbolic limit cycle 
with ri = r and 6i = 6j for all i and j. If the limit cycle is stable the truncated 
system is synchronized. Eigenvalues (38-39) are expressed in terms of the coupling 
load parameters, and implicitly define the region corresponding to stable synchronous 
behavior in parameter space. 

It remains to show that the stability of the synchronous solution in the truncated 
system implies the existence and stability of nearby synchronous solution in the 
original system (15-16). Let Aj_ij = 6j — 6j_i for j = 2,...,N, and let x = 
(Ai^2, • • • 7 '^N-i,N, ri, ■ . . , ttv). The only 0{1) terms in equations (34-35), are the unit 
terms in (35). By definition of Aj_ij, these terms do not occur in the differential 
equation for x- K follows that that in the new coordinates {x,Oi), equations (34-35) 
have the form 

X' = eF,{x) + e'F2{x,0i) 

e[ = l + eG,{x) + e'G2ix,0,). (40) 
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The explicit form of Fi, F2, Gi, and G2 is not of importance. 

To study the stabihty of the synchronous state Xo = (0)'')) we note that the 
eigenvalues of eD^Fi{xo) are given by (38-39) when r = (r, r, . . . , r) and r satisfies 
(36). The following proposition shows that if the eigenvalues Aj have non-zero real part 
then they completely determine the stability of the synchronous state for the full system 
(15-16). 

Proposition 2 If Fi{xo) = 0, and D^Fi{xo) has eigenvalues with non-zero real part in 
(40), then there exists an eo, such that system (40) has a limit cycle 0{e) close to the 
limit cycle of the unperturbed system, obtained from (40) by setting F2 = G2 = 0, for 
all e < eo . The Floquet exponents of the perturbed limit cycle agree with the eigenvalues 
ofD^Fixo) to Oie). 

Proof: This is a consequence of standard results about near identity changes of 
coordinates for systems with a single frequency. See [21, Chapter 1.22] and references 
therein. □ 

Proposition 3 Assume that a nonzero solution r of equation (36) exists. If the 
eigenvalues given in (38-39) have negative real part, then (15) has a stable, synchronous 
periodic solution. If one of the eigenvalues has positive real part, this periodic solution 
is unstable. 

Proof: As a consequence of Proposition 2, depending on the real part of the 
eigenvalues in (38-39) there is either a stable or unstable synchronous solution of system 
(25-28). Since the change coordinates Bw = z affects all of the pairs of coordinates 
{wkiWk) in the same way, a synchronous solution in the w coordinates corresponds to 
a synchronous solution in the original z coordinates. The stability properties of this 
solution are preserved under a linear change of coordinates. □ 

Remark 5 Theorem 1 only states that the truncated normal form provides an 0{e) 
approximation on timescale of 0{l/e), and typically this approximation does not hold 
on longer timescales. Nevertheless, Proposition 2 implies that the limit cycle of the 
truncated normal form approximates the limit cycle of the original system to 0{e), since 
the approximation of the amplitude x is valid for all time. 

Remark 6 In ['L">], a stroboscopically discretized system was used in a heuristic 
argument to obtain similar result in the case of van der Pol oscillators. There it was 
necessary to determine the angular frequency of the inphase solution to 0{e) in order 
to find appropriate strobing period. The angular frequency of the periodic solution can 
also be estimated from (34-35) up to second order as 

^ = l-e(e(r) + ^sin5). (41) 
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Figure 2. Floquct multipliers for a van der Pol oscillator array. Solid line represents 
analytical result, while dots represent results obtained from numerical calculations. 
The coupling parameters are k = 1, e ~ 0.1, and (a) fi2 ~ 0.8, (b) fi2 = 1-5. 



5. Examples 

The synchronization condition for the array of globally coupled oscillators is obtained 
from (38-39) by setting Aj < 1 for alH > 1, and can be expressed in terms of control 
parameters /ii and fi2- In general, our method allows us to do a single calculation 
for a certain network configuration, and then obtain results for a variety of different 
oscillators in a straightforward fashion. To find the synchronization condition for a 
specific oscillator type it suffices to find the function (p^ which characterizes the resonant 
terms in the equation of motion of the uncoupled oscillator. We illustrate this in several 
examples and compare our approximating solutions with numerical results. Agreement 
between analytical and numerical calculation is good even when using only first order 
corrections in e. We retrieve result from [15] as a special case of our general result. 



5.1. Van der Pol oscillator arrays 

Consider a network of globally coupled van der Pol oscillators described in Sec. 2.1. 
From (14) we find that the resonant terms are described hj (p^ = 1/2 — wiv/8, and 
hence R{r) = 1/2 — 6(r) = 0. From (36) we find that the inphase state exists for 
r = 2^1 + Kcos6/Z. By substituting expressions for R(r) and 9(r) in (38) and (39) 
we find the Floquet exponents for the synchronous solution 

A2 =-|(l + |^o«^) (42) 

An,n+i = - ^ (l + 2| COS 5) ± |y(l + |cos5)'- (|sin5)'. (43) 

In order to test our results we evaluate Floquet multipliers || numerically and compare 
them to the values estimated from (43). The results are 0{e) close (Fig. 2). 

II In these comparisons we use Floquet multipliers for convenience, since they are easier to handle 
numerically. The Floquet multipliers are given by = ^^'^^d^ _|_ C(e^), where 10 is given by (41). 




Figure 3. Floquct multipliers for a van der Pol-Duffing array. Coupling parameters 
are a = —0.2, k = 1, e = 0.1, and (a) \i2 = 0.8, (b) \xi = 1.5. 



5.2. Van der Pol-Duffing equation 

In a recent study of micromechanical and nanomechanical resonators models using 
parametrically driven Duffing oscillators [2")] and van der Pol-Duffing oscillators [4] 
are proposed. The van der Pol-Duffing equation 

x" + x- e(l - - eax^ = 0, (44) 

is obtained from the van der Pol equation by an addition of a cubic term. By switching 
to complex coordinates (11) and writing (44) in normal form we obtain 

z' = — iz + - {4:Z — z^z + iZaz^z) 
8 

z' = iz — - {—Az + zz^ + i'iazz^). 
8 

The resonant part of the nonlinearity is given by ^{f) = (4 — + i3ar^)/8. The real 
part of (f)^ is the same as in the case of van der Pol oscillator, so this system has the 
same inphase solution. Substitute the real part R{r) = 1/2 — r^/8 and imaginary part 
6(r) = 3ar^/8 of (p^ in (38) and (39) to find the approximate Floquet exponents 

As = - ifl + ^cos5) (45) 



2\ Z 

A„,„+i = - |(l + 2|cos5) (46) 



±-Ml + -cos6j -(^-sin^j + 3a-(2sin5+ -sin2(5). 

The numerical simulations (Fig. 3) support our result. Note that it follows from (20) 
that our approximations can be expected to break down for small values of /xi. 



5.3. Synchronizing sources and sinks 



As noted in Remark 4, it is possible to turn a network of systems with flows that have 
only a source (or sink) at the origin into a network of synchronous limit cycle oscillators 
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Figure 4. The radius of the periodic, inphase sohition is obtained by solving 
(36). For an array of elements (47) it is given by the intersection of the curve 
rR{r) — 7'/2(r^ — 3.6r + 3.32) and the line — k cos (5/(2Z)r. If the slope of the line 
is too small only the trivial solution r = exists, and the system does not exhibit limit 
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Figure 5. Phaseplane (left) and time diagram (right) of inphase limit cycle solution 
for coupled sources (47). Dashed line represents the solution for the truncated normal 
form system, and solid line numerical solution for the full system. 



with an appropriate coupling. Consider a dynamical system described by 

3.6 



X 



X — ex 



Ax' 



X 



3.32 



0. 



(47) 



Each individual system has only an unstable fixed point at the origin, and no other 
repellers or attractors. Coupling these systems as in (15-16), and switching to complex 
variables (11), gives the following equations of motion 

z' = -iz + ief{z,z) (48) 
z' = iz-ief{z,z) (49) 

where the nonlinear term is 



f{z,z) 



2i 



(z + z) 



3.6 



3.32 



(50) 
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Figure 6. Floquct multipliers for coupled sources, k — 1, e ~ 0.1, and (a) fi2 ~ 0.8, 
(b) fi2 - 1.1. 

If we keep resonant terms only, (48) and (49) become 

z' = -iz + ez(j)^{zz) (51) 
z' = iz + ez(j)^{zz) (52) 

with (p^ = R{r) = l/2{zz - 3.6v^ + 3.32). 

From (36) we find that the hmit cycle solution does not exist for —kcos6/{2Z) < 
0.08. Below that value the coupled elements behave as unstable foci, which can 
be easily checked numerically. As —kcos6/{2Z) increases above 0.08 the system 
undergoes a supercritical saddle-node bifurcation of limit cycles, in which a stable and 
an unstable inphase limit cycle are created. These two solutions are are represented as 
the intersection of the curve R = r/2(r^ — 3.6r + 3.32) and the line R = —K,cos6/{2Z)r 
in Fig. 4. For a suitable choice of coupling parameters it is possible to obtain inphase 
limit cycle solutions for the array of sources. In Fig. 5 we show oscillations of an element 
in the array, when the coupling parameters are set to k = 1, /ii = 1.2 and fi2 = 0.8. The 
Floquet exponents for both limit cycles are readily obtained from (39). These results 
agree with numerical calculations to the expected error. In Fig. 6 we present results for 
the "stable" cycle. 

We note that a direct extension of the method discussed in the previous sections 
was used to derive these results, since the nonlinear term in (47) is not a polynomial 
(see [22]). 

6. Conclusion 

In our study of synchrony in globally coupled oscillator arrays we introduce a normal 
form based method, which has a number of advantages over methods commonly found in 
the physics literature. The method provides a clear and mathematically rigorous way of 
finding the onset of synchronization in the array with respect to changes in the control 
parameters. It allows for an easy distinction between contributions to the dynamics 
coming from the network configuration and those coming from the internal structure 
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of the network elements. As a consequence, we are able to carry out calculations 
for particular network architectures without having to specify the exact form of the 
nonlinearities in the individual elements. 

To apply these results to specific weakly nonlinear oscillators it is only necessary 
to find the function (j)^, which characterizes the resonant part of the nonlinearity, and 
substitute it in the general solution thus obtained. The function (j)^ does not depend on 
the coupling scheme, and is easily derived from the equations for the uncoupled systems. 
It is therefore tempting to think of synchronization classes as different systems may lead 
to the same (f)^. 

Although we have chosen a very specific linear coupling in our exposition, the 
method can be applied to a variety of network configurations simply by retracing the 
steps we outlined. A similar calculation can be carried out even if there is a weak 
nonlinearity in the coupling equation (16) itself. Moreover, the method can be extended 
to higher orders in the small parameter in a straightforward fashion, and the procedure 
can be automated. 

The analysis of the synchronous solution of the network (15-16) was particularly 
simple due its Sn symmetry. When the network has less symmetry, or only local 
symmetries (see [26]), a similar reduction can be performed. In such networks one 
expects polysynchronous solutions, in which groups of oscillators within the network 
oscillate synchronously. One can further expect to obtain a pair of equations for each 
cluster of oscillators in the network [12]. The stability of these clusters can then be 
analyzed in a manner similar to the one introduced in this paper. 

Although we have not treated the case of Josephson junctions, a similar analysis 
can also be performed. In fact we believe that coherent behavior in networks of various 
dynamical systems can be studied using the method of normal forms, and intend to 
investigate this in the future. 
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